Propagation of Initially Excited States in Time-Dependent Density Functional Theory 
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Many recent applications of time-dependent density functional theory begin in an initially excited 
state, and propagate it using an adiabatic approximation for the exchange-correlation potential. This 
however inserts the excited-state density into a ground-state approximation. By studying a series 
of model calculations, we highlight the relevance of initial-state dependence of the exact functional 
when starting in an excited state, and explore the errors inherent in the adiabatic approximation that 
neglect this dependence. 



I. INTRODUCTION 

Time-dependent density functional theory (TDDFT) 
is an exact reformulation of the time-dependent quan- 
tum mechanics of many-electron systems [1, 2] that 
operates by mapping a system of interacting elec- 
trons into one of fictitious non-interacting fermions, the 
Kohn-Sham (KS) system, reproducing the same time- 
dependent density of the true system. By the theorems 
of TDDFT, all properties of the true interacting system 
may be obtained from the KS system, in principle. In 
practice, approximations are needed for the exchange- 
correlation effects: in particular, for the exchange- 
correlation (xc) potential, Vxc[n'-'i'{0),^{0)]{r,t), a cru- 
cial term in the single-particle KS equations determining 
the evolution of the ficitious non-interacting fermions. 
This functionally depends on the one-body density, 
n(r, t), the initial many-body state of the interacting sys- 
tem \I/(0), and the initial state of the KS system $(0). 

From the inception of the theory, subtleties were 
noted in the functional-dependences, which make good 
approximations for the xc potential more challenging 
to derive than for its counterpart in the (much older) 
ground-state density functional theory [3-5]. One of 
these is memory: the dependence of the potential at time 
t on the density at earlier times, and on the initial-states 
\l/(0) and $(0). Efforts to include memory-dependence 
in functionals [6-14], are not widely used, and have 
all focussed on the dependence on the history of the 
density, neglecting the dependence on the initial states. 
The vast majority of applications have been in linear re- 
sponse from a non-degenerate ground-state, and there 
this initial-state dependence (ISD) is redundant: by the 
Hohenberg-Kohn theorem, a non-degenerate ground- 
state is a functional of its own density. In fact, the simple 
adiabatic approximation, where the instantaneous den- 
sity is input into a ground-state approximation, neglect- 
ing any memory-dependence, has shown remarkable 
success for a great range of excitations and spectra. All 
commonly available codes with TDDFT capabilities are 
written assuming the adiabatic approximation. Over the 
years it has been understood that certain types of exci- 
tations cannot be captured by the adiabatic approxima- 
tion, e.g. double excitations [15], excitonic series in op- 
tical response of semiconductors [16], and there is both 



on-going research in developing frequency-dependent 
kernels to deal with this, as well as understanding when 
not to trust the calculated adiabatic TDDFT spectra. 

With the results and understanding of recent years 
lending a level of comfort with calculations of spec- 
tra and response, TDDFT has now entered a more ma- 
ture stage, and with that, comes more adventurous ap- 
plications. In particular, for electron dynamics in real 
time, evolving under strong laser fields, or coupled 
electron-nuclear dynamics following photo-excitation 
(e.g. Refs. [17-19]). The role of memory is less well- 
understood in these applications; studies on model sys- 
tems have shown that sometimes memory-dependence 
is essential [20-24], other times it is not important at 
all [25, 26]. Moreover a new element enters: the initial- 
state dependence that was conveniently and correctly 
brushed aside in the linear response regime, now raises 
its head. In applications such as modeling solar cell pro- 
cesses, the initial photo-excitation of the electronic sys- 
tem is not dynamically modeled; instead the dynamics 
begins with the electronic system assumed to be in an 
excited state. The initial state is not the ground-state, yet 
no truly initial-state dependent functionals are available 
today, hence we ask how large are the errors in such cal- 
culations? Knowing that the true interacting system be- 
gins in a certain excited state, is there an optimal choice 
for the initial KS state when propagating with an adia- 
batic approximation? How large are the errors due to 
ISD compared to those due to history-dependence? 

In this paper, we begin to answer these questions 
by considering a series of model calculations of two- 
electron systems. We start by reviewing the underlying 
theorems of TDDFT and the subtleties of ISD, even in sit- 
uations where we may not expect it. We show how ISD 
leads to a non-zero xc potential even for non-interacting 
electrons. 

Then, we move to electron dynamics in the 
model soft-Coulomb helium atom, performing adia- 
batic TDDFT calculations with different initial-states 
and comparing with exact dynamics. Here we must 
dissect other sources of error in usual TDDFT calcula- 
tions, such as using an initial KS excited state whose 
density does not equal the true excited state density, in 
order to properly ascribe the influence of ISD. Finally 
we make the first step towards investigating how ISD 



can affect an area currently of much interest, namely 
coupled electron-ion dynamics, by performing Ehren- 
fest calculations for a model LiH system. We work in 
atomic units throughout (e^ ~ nie = h ^ 1). 



orbitals (j)i is selected, with the required property that 

N 

J2\Mr,0)\' = n(r,0)aiid (3) 

N 

-V-Im^0:(r,O)V(/),;(r,O) = n(r,0) (4) 



where 



II. INITIAL-STATE DEPENDENCE IN TDDFT 



In TDDFT, one evolves a set of single-particle orbitals 
{(?!)j(r, t)} with a one-body KS potential: 



d 



--V' + Vs{r,t)]<P^{r,t) 



(1) 



vs{r,t) = ^;ext(r,i)+WxcM(r,i)+Vxc[?^;*o,*o](I■,i) (2) 

wherewH(r,t) — /ri(r', t)/|r— r'jd'^?'' is theusualHartree 
potential and Uxc(r, t) is the xc potential which depends 
on the entire history of the density, the interacting initial 
state 4*0, and the initial KS wavef unction $0. 

The origin of the initial-state dependence in the func- 
tionals is that the one-to-one Runge-Gross mapping [1] 
between time-dependent densities and potentials holds 
for a fixed initial-state. Thus, the external potential Vc^t 
functionally depends on the density and the true ini- 
tial state '5(0), and the KS potential v^ depends on the 
density and the KS initial state $(0). These functional 
dependences are not directly relevant themselves in a 
practical calculation, because there only the xc potential 
Vxc needs to be approximated. However, they lend their 
dependences to Wxc via Eq. (2), which must therefore de- 
pend on both initial states and the density. 

There are two known situations in which there is 
no initial-state dependence. When the initial state is a 
ground-state, then, by the Hohenberg-Kohn theorem of 
ground-state DFT, it is itself a functional of its own den- 
sity. ISD is redundant, as the information about the ini- 
tial state is contained in the initial density. (See also 
Sec. II A). The other situation is for one-electron sys- 
tems: starting in any initial state, there is only one poten- 
tial that can yield a given density-evolution [27]. In all 
other situations, it is assumed that ISD cannot be sub- 
sumed into a density-dependence; explicit demonstra- 
tions for two electrons can be found in Refs. [22, 27]. 
Ref. [28] derived an exact condition relating the depen- 
dence on the history of the density to ISD. This condition 
is likely violated by any history-dependent functional 
approximation that has no ISD. 

Technically, one may choose any initial KS wave- 
function that has the same n and ii as the true initial 
state [29]. Usually a single Slater-determinant of A^ spin- 



n{r,0) = nY^ I dx2--- Ma;Af|*(a;,X2 • • • xjv, 0)^ (5) 
ri(r,0)= -iVV-Im < y^ dx2 ■ ■ ■ / dxN'^*{x,X2XN,0) 

\/^{x,X2,,XN,0)> (6) 

where x = (r, a) and J dx = J^a I '^^^ ■ 

In this paper we will consider different choices of $ (0) 
for two-electron systems that begin in the first excited 
singlet state of the true system, denoted 4'*. Specifi- 
cally we will, at various points, investigate the follow- 
ing three forms: 

(a) An excited KS singlet state with two occupied or- 
bitals 00 and 01, whose spatial part has the form 

$*(ri,r2) = -^(</)o(ri)0i(r2)+0o(r2)0i(ri)) . (7) 

Note that this state is a sum of two Slater determinants. 

(b) A ground KS singlet state, $^'', which, for two elec- 
trons corresponds to a single doubly-occupied orbital 



<i>s^K](ri,r2) = (/.(ri)0(r2), 



with 



r) = ^Jn*{r)/2 



(8) 



(9) 



where n*{v) = 2/ \^*{v,v')\^d^r' is the density of the 
initial excited state. This is a single Slater-determinant, 
and is the usual choice when the true initial state is a 
ground state, 
(c) A spin-symmetry-broken excited KS state of the form 



*SB(i-V,ra) 



1 

V2 



/)o(r)(5o-t 0o(i"')'5o-'t 

hho)5ai (l)l{v')5a'i 



(10) 



This state is not a spin-eigenstate (although it does have 
(S) = 0), but is a valid choice for an initial KS state, hav- 
ing the same total density as the true system. The up and 
down spin-densities are not equal, and we shall evolve 
them in different spin-up and spin-down KS potentials, 
but, unlike in spin-DFT, we do not consider the spin- 



densities separately meaningful; only their sum is con- 
sidered as an observable [30]. Instead, (c) will be used 
to illustrate the effects of orbital-specific functionals, as 
will be discussed more later. 

Now the exact KS potential differs depending on 
which choice (a), (b), or (c), is made, but any adiabatic 
approximation is identical for them all. Almost all the 
calculations being run today, certainly all that are coded 
in the commonly available codes, utilize such an ap- 
proximation, and insert the instantaneous density into 
a ground-state approximation: 



<-[n;vI/(0),$(0)](r,i)=<c[«W](r) 



(11) 



It is not surprising that there are errors inherent in such 
an approximation but one question we hope to shed 
light on by our investigations here, is how much of this 
error is due to the lack of ISD, rather than the lack of 
history-dependence (dependence on n{r,t' < t)). Only 
the latter occurs in usual calculations that begin in an 
initial ground-state, where, as mentioned earlier, ISD 
may be subsumed into density-dependence. 



A. Beginning in the ground-state: ISD or not? 
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FIG. 1: Top Panels: The first excited-state density, n*{x), 
for two non-interacting electrons in a one-dimensional har- 
monic potential (left) and soft-Coulomb potential (right). 
Lower Panels: The corresponding correlation potentials, 
uc[n*, **, $^''[n*]](a;), for an initial KS wavefunction chosen 
to have ground-state form but yielding the initial excited-state 
density n* (case (b) of text). The insert shows the KS potential 
(solid line) compared with the external potential (dashed line), 
chosen as harmonic and soft-Coulomb, respectively. When the 
initial KS state is chosen to be <I>* (case (a) in text), on the other 
hand, the correlation potential is zero. 



Many practical situations start with the system in its 
ground-state; in fact most calculations, except in the 
most recent years, have assumed an initial ground-state. 
This is fortunate for TDDFT, since, at least at short times, 
the adiabatic approximation with its ground-state func- 
tionals, which have become increasingly sophisticated 
over the years, could then be expected to work reason- 
ably well. However, we will show that the the subtleties 
of initial-state dependence can appear even in this case. 

A completely legitimate choice for the KS initial state 
would be the true ground-state wavefunction as it triv- 
ially satisfies the initial conditions Eq. (3) and Eq. (4). 
The exact KS potential, for this choice, reproduces the 
density evolution of the exact wavefunction propagat- 
ing in the interacting system by propagating the exact 
wavefunction in a non-interacting system. 

Given that at the initial time we are dealing with a 
ground-state, one might expect that the adiabatic ap- 
proximation would be exact, at least initially, however 
this is not the case. The adiabatic approximation re- 
turns the xc potential for a ground-state DFT calcula- 
tion, where the KS wavefunction is the ground-state KS 
wavefunction. If we were to start the KS calculation 
with this ground-state KS wavefunction, the adiabatic 
approximation is then exact at the initial time. 

However the initial interacting wavefunction is not a 
ground state of the non-interacting KS system, and so, 
choosing this as the initial KS state means the adiabatic 
approximation will be in error from the start. 

Thus care must be taken when using the short-hand 
that there is no initial-state dependence when starting in 
the ground-state. More precisely, what is meant is that 



1) the initial KS orbitals are chosen to be those of the 
KS ground-state wavefunction, and 2) the dependence 
on initial state for this case can be subsumed into the 
density via the theorems of ground-state DFT. 



III. NON-INTERACTING ELECTRONS 

The importance of initial-state dependence is strik- 
ingly evident even for the hypothetical case of non- 
interacting electrons. Due to ISD, the xc potential is not 
always equal to zero, which may be surprising at first 
sight, given that there is no interaction. Although no- 
one in their right mind would perform KS calculations 
for non-interacting electrons for practical purposes, it is 
instructive to consider how the KS system behaves in 
this case. In particular, the studies suggest what is the 
best choice of KS initial state for a given true initial state 
when an adiabatic approximation is used. Since the adi- 
abatic approximation is designed for ground-states, is 
the error least if we always choose a KS state that is a 
ground-state? 

In the following we consider the "true" system to be 
non-interacting, i.e. we scale the electron-electron inter- 
action by A in the limit that A ^- and consider terms 
zero-th order in A only, e.g. the Hartree potential van- 
ishes. We then consider the xc potential when the KS 
system is started in different allowed initial-states. 

Consider two non-interacting electrons prepared in 
an excited state ^* and evolving in some potential 
VcKt(t)- To start the KS evolution, any initial state of the 



same density n, and first time-derivative, n, zero in this 
case, may be chosen. We consider the initial state choices 
(a) and (b) introduced in Sec. II, which become here: (a) 
$(0) = **, and (b) $(0) = Mr)Mr') = *^% <^o = 
■\/n*/2, where n* is the density of ^*. 

For choice (a), the exact xc potential vanishes, 
Wxci'T-; ^*, *i'*](r, t) = 0. This can be most easily seen by 
invoking the uniqueness property of the Runge-Gross 
mapping: for a fixed particle-interaction, only one po- 
tential {vcxt{t)) can yield a given density-evolution from 
a given initial state. Since $* = ^* for choice (a), we 
conclude Uxc(0 ~ 0- This result can also be seen from 
the general formula [29]: 

V-[n(r,i)VK(r,i)+z;xc(r,i))]=9(r,i)-'7s(r,t) (12) 
where 

q(r,i) =lv.(vl/(t)|[j(r),f + Kc]|*W> (13) 
in 

qs{r,t) =lv.($(i)|[j(r),T]|$(i)) (14) 

and j (r) is the current-density operator. When the two 
initial wavefunctions are the same, the RHS vanishes at 
i = in the non-interacting limit, and we are left with 
Vnit = 0) = 0, and Vxc{t = 0) = 0. Stepping forward in 
time, it follows that Wxc(0 ~ for all times. 

Now consider choice (b). In this case, 

v^c[n; **, $ss](r, 0) is non-zero: the RHS of Eq. (12) at 
i = is no longer zero in the non-interacting limit, due 
to the difference in the initial wavef unction in q and Qs . 
That this is non-zero is also expected by the fact that v^c 
must be such that Vs = Woxt + Vhxc evolves the initial 
$s^ with the same density for all time as the different 
initial state ^* has when evolved in Vcxt- So, ISD leads 
to a non-zero xc potential, even when the electrons do 
not interact. 

For an explicit demonstration, consider the initial 
time and realize that Wxc(i = 0) depends entirely on the 
initial states, as can be seen from Eq. (12), and not on 
the choice of the external potential. Figure 1 plots this 
potential v^dt ~ 0) for the first excited state of an ex- 
ternal harmonic potential {^x^, on the left), and a soft- 
Coulomb potential {—2/\/jp~{~T, on the right) in one- 
dimension. The top panels show the density, and the in- 
sets show the external potential in which 'J* is an eigen- 
state, as well as the the KS potential Wcxt + "c in which 
^^^ is an eigenstate (the ground-state, since there are no 
nodes). Because for two electrons in a spin-singlet with 
one doubly-occupied orbital, fx = — Wh/2 = in the 
non-interacting limit, this effect appears entirely in the 
correlation potential, and can be interpreted as static cor- 
relation: due completely to the non-single-determinantal 
structure of the true initial state. In case (a), this effect 
vanishes, because the KS initial state is chosen also not 
to be a SSD. 

We now ask what an adiabatic potential would give: 



approximating Vxc[n; *(0), $(0)] by a ground-state po- 
tential w|c['^]- To distinguish errors arising from the 
choice of approximate ground-state functional itself, 
we consider an "adiabatically exact" potential [21, 25], 
'^xc'^ ''^["l- ^6 shall define this generally in the next sec- 
tion, but for now it suffices to define it such that if both 
the true and KS wavefunctions at all times were in fact 
ground states of some potential, then the exact xc po- 
tential is the exact ground-state one, Vxc[n; "if^^, ^^^] ~ 
v^^^^^™[n]. Consider again just the initial time, t = 
0. In the non-interacting limit, for both the true and 
KS wavefunctions to be ground-states that have the 
density of the excited state **, then 'I'(O) = $(0) = 
^no(r)no(r')/2, and we rapidly conclude 
that 

^;xc[n,*^^$^^](r,i) = vfj^-''^[n]{r,t) = (15) 

in the non-interacting limit. Applying now this adia- 
batic approximation to the case of the initial true excited 
state ^* considered above, we conclude that the adiabatic 
approximation is exact for choice (a), when the initial KS 
state is also chosen excited, as we had argued there the 
exact Vxc also vanishes. It is inaccurate for choice (b), 
when the initial KS state is chosen to be a ground-state 
one, where the exact correlation potential is non-zero. 

This study suggests that in the general interacting 
case, errors in adiabatic TDDFT will be least when an 
initial KS state that most closely resembles the configu- 
ration of the true excited state is chosen. This expecta- 
tion is indeed borne out in the following studies. 



IV. MODEL TWO-ELECTRON SOFT-COULOMB 
INTERACTING SYSTEMS 

The case of non-interacting electrons illustrated the 
importance of ISD, while offering a hopeful diagnosis 
for the adiabatic approximation; the latter becomes ex- 
act if the initial state is chosen appropriately. Of course, 
electrons do interact, and now we turn to studying the 
effect of ISD on the dynamics of interacting systems. 

When running an approximate TDKS calculation 
starting in an excited state of the interacting problem, 
three separate sources of error come into play. First, 
excited states of the exact ground-state KS potential 
tis[7io](r) do not have the same densities as interacting 
excited states of the potential Wcxt (r) whose ground-state 
density is no(r). Yet these are usually the ones chosen in 
practice. We discuss this problem in Section IV A. The 
second source of error is the central one for this paper: 
the use of the adiabatic approximation, when, from the 
start, we have an excited state. We consider first the 
"adiabatically exact" approximation, mentioned also in 
Sec. Ill, but considered now for interacting electrons, in 
Section IV B. This will allow us to separate the errors 
from the choice of ground-state functional approxima- 
tion used in the adiabatic approximation, which is the 




FIG. 2: Exact excited-state density (solid line) compared to the 
exact excited-state density of the exact ground-state KS poten- 
tial (dashed line) for the soft-Coulomb Helium atom. 



third issue. In Sections IV D, and V, we will study the 
effect that missing initial-state dependence has in prac- 
tice, on electron dynamics in model two-electron sys- 
tems when begun in the first singlet excited state of the 
system. 

In our model one-dimensional Helium atom, the two 
electrons interact via a soft-Coulomb electron-electron 
interaction, Vc 



an external soft-Coulomb potential. 



v^^tix) = -2/Vx^ + l 



1 and live in 



(16) 



Such a model, straightforward to solve numerically, is 
popular in understanding and analyzing electron inter- 
actions, in both the ground-state and in strong-field dy- 
namics [31], inside and outside the density-functional 
community. In our model of the diatomic LiH molecule, 
the soft-Coulomb interaction is also used, while the ex- 
ternal potential has an asymmetric soft-Coulomb well 
structure (Sec. V). 

Before proceding, we make a computational note. The 
true dynamics are propagated on a real-space grid using 
the exponent midpoint rule Taylor-expanded to fourth 
order. Imaginary time propagation along with Gram- 
Schmidt orthogonalisation is first performed to find 
the lowest eigenstates. The adiabatic exact-exchange 
(AEXX) dynamics uses the Crank-Nicolson method 
with an explicit predictor step for the Hartree-exchange 
potential. In both cases, a timestep of 0.001 is used with 
a grid spacing of 0.1. When possible these calculations 
were tested for accuracy against the parallelized code 
OCTOPUS[32], which was also used for LDA and LDA- 
SIC runs. 



A. The Excited-State Density 

The theorems of TDDFT require that the KS initial 
state has the same 7i(r,0) and 7i(r,0) as the true inter- 
acting system (Sec. II). In this section we discuss the 
difficulty of fulfilling this requirement when the initial 
state is not a ground-state. 




FIG. 3: Top Panel: The exact xc potentials for the two choices 
of initial state, "3?* (solid line) and $^^ (dot-dashed line). Lower 
Panel: The KS potentials for which the two initial state choices, 
<&* (solid line) and <1>^^ (dot-dashed line) are eigenstates. Also 
shown are the exact ground-state KS potential (dashed line) 
for the soft-Coulomb Helium potential (dotted line). 



When the true initial state is a ground-state, the nat- 
ural and usual choice for the KS initial state is the non- 
interacting ground-state, and this can be found by solv- 
ing the ground-state KS equations. However, if we want 
to compute the dynamics of an excited state, we en- 
counter the problem that there is no DFT scheme to find 
excited state densities[34]. Furthermore, even if we use 
a more computationally expensive higher-level wave- 
function method to calculate the state 'I'(O), we still must 
then choose an initial $(0) in which to start the KS cal- 
culation. If we start the interacting system in a first ex- 
cited state of some t'cxt(i"), the results of the earlier sec- 
tions suggest a good choice for adiabatic TDDFT is to 
start the KS system in a corresponding non-interacting 
excited state where one electron is excited from the high- 
est occupied to lowest unoccupied orbital of some Vs 
(as in form (a) of Sec. II). This 'Us(r) however cannot 
be the ground-state KS potential fs(r) corresponding 
to the interacting ?;oxt/ because this would not satisfy 
Eq. (3): Vs{r) yields the same ground-state density in a 
non-interacting system that Wcxt (r) yields in the interact- 
ing one, but the densities of their excited-states are dif- 
ferent. Yet, the usual practice is to use the correspond- 
ing excited-state of the ground-state KS potential [33]. 
This is partly because the density of the interacting ex- 
cited state is not often known anyway In Fig. 2 we 
compare the exact density of the first excited interact- 
ing singlet state in the soft-Coloumb He atom (Eq. 16) 
with the density of the excited state of the correspond- 
ing ground-state KS potential. Although matching the 
general shape, the latter is too narrow and misses some 
structure. 



To be able to separate the error from not quite hav- 
ing the exact initial density from the error from using 
an adiabatic approximation, we now search for a non- 
interacting system where the first singlet excited-state 
density is exactly equal to the true excited-state den- 
sity n*, for our model He atom. Several such potentials 
may exist [34], and in the lower panel of Fig. 3, we plot 
one such KS potential, denoted Vs[n* , ^*]{x) and com- 
pare to the KS potential where n{x) is the exact ground- 
state density v^ [n] {x, 0) of the model He atom. (The first 
excited state densities of these potentials are precisely 
those shown in Fig. 2). We also show Vs[n*]{x,0), i.e. 
the KS potential for which n*{x) is the density of the 
ground-state. This corresponds to choice (b) in Sec. 11 
for the initial KS state. 

Once a valid KS state is found satisfying Eqs. (3) 
and (4), the initial xc potential is completely determined: 
Eqs. (12) and (13) show that at the initial time, v^c is a 
functional of just the two initial states ^(0) and <i>(0). 
This can be added to any external potential, prescribed 
by the physical problem at hand, and Hartree-potential, 
determined by the initial density, to start the evolution. 
In other words, the choice of the initial KS state funda- 
mentally points to an xc potential, rather than to a KS 
potential, and so in the top panel of Fig. 3 we plot the 
xc potentials corresponding to choices (a) and (b) of the 
initial KS state. In the following examples, we take the 
initial v^xt to be that of the soft-Coulomb Eq. (16); some- 
times subsequently an external field is added. 

We have therefore, now found initial states of the form 
(a) and (b) of Sec II, both yielding the same density 
as that of the true excited-state, meeting the conditions 
Eq. (3) and (4). We can now move forward to investigate 
the time-dependent properties. 



B. The Exact Adiabatic Approximation 



v^(,[ii ,V ,* ](x,0) 
Vj,j,lii*.>P*,*^" J(x,0) 
'"' " [n*l(x,0) 




Vj,[n ,<!> ](x,0) 
Vj,[n*,<I>'^'](x,0) 
'■"■' (x,0) 



FIG. 4: Top PanebThe exact xc potentials and the exact adia- 
batic xc potential. Lower Panel: The exact and adiabatically- 
exact KS potentials. 
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We next move to the error introduced by using the 
adiabatic approximation. 

To isolate this error we will calculate the adiabatically 
exact potential at the initial time, and start in a KS state 
of the exact same density as the true state, as found in 
Section IV A. 

The adiabatically exact potential is defined by [25] 



«5^'^-'='^W 



,,adia— all r 



, aclia— all r 



n - f H n 



(17) 



where 



,,adia— all 



[n](r,<) is the potential for which 



FIG. 5: The density at times 5au, lOau, 15au, and 20au for ini- 
tial states $* (solid line) and $^^ (dashed line) propagating in 
the static exact-adiabatic initial potential. Also shown is the 
exact excited-state density (dotted line) which is the exact so- 
lution at each time. 



doubly-occupied orbital, 0(r) = ^n*(r)/2 where n*{r) 
is the density of the initial state 4**. 



The interacting potential, v, 
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i](r,0) is solved for 



n(r,t) is the non-interacting ground state density and 
v^^l'^~^^^[n\{x) is the potential for which interacting elec- 
trons have ?7(r, t) as their ground-state density. If both 
the true and KS wavefunctions are actually always in 
some ground-state, then Eq. (17) becomes the exact xc 
potential. We shall consider this potential for initially 
excited states, only at the initial time. 

The non-interacting potential i'g'^'^~''"[n*](r, 0) may 
be easily found by inverting the KS equation for a target density, we have found u, 



using an iterative technique whereby the external po- 
tential is updated based on the difference between the 
ground-state density of the current iteration and the 
density we are targeting. This is based on the inversion 
algorithm of Ref . 35, but generalized to the interacting 
case[25]. We also increase the update to the potential in 
regions of low density by simply using the inverse of the 
density (up to a maximum value) as a weighting factor, 
similar to Ref. 30. When the density converges to the 
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n\{x). 



The dotted line in the lower panel of Fig. 4, shows 
at the initial time, the full KS potential with the adia- 
batically exact xc potential, namely v^'^^^{x) = Uoxt(a;) + 
VH[n]{x) + v^'^^^^^^[n]{x), choosing Vcxt{t = 0) as the soft- 
Coulomb potential of Eq. 16, for the first excited singlet 
state ^*. Alongside, we compare this to the exact KS po- 
tentials found in Sec. IV A (see Fig. 3), for the initial state 
choices of case (a) and (b) respectively. The adiabatically 
exact KS potential ^^'^''^(x) tracks the shape of the ex- 
act KS potential Us [n* , ^* , $*] although misses structure. 
The difference between them is the difference in their xc 
potentials shown in the upper panel of Fig. 4. The dif- 
ference is much larger for the comparison with the exact 
KS potential when the initial KS state is chosen as the 
ground-state, (i.e. v^'^^^~°^[n*] is closer to v^c[n*, ^*, $*] 
than it is to Vxcin*,"^*, ^^^[ii*]]); in this case, there is a 
large static-correlation contribution to the exact xc po- 
tential. The results here are consistent with the non- 
interacting case discussed in Sec. Ill: with interaction, 
the adiabatic approximation is no longer exact for the 
choice of $*, however it is still the better choice over 

The fact that the adiabatically exact potential does not 
equal the true potential, even at the initial time, leads to 
erroneous dynamics. To illustrate this effect we propa- 
gate our exact initial wavefunctions in the adiabatically 
exact potential and compare with exact dynamics. Since 
calculating at each time is computationally demanding 
and delicate, we will simply evolve in time without any 
additional perturbing potentials, and hold the adiabatic 
potential fixed to its initial value. The idea is that since 
the exact evolution is static, the adiabatically exact po- 
tential is static and equals its initial value at all times; if 
this was a good approximation, it would yield density- 
dynamics that were close to static. The deviation of the 
adiabatic propagation from the exact is a measure of its 
error. Note such a calculation does not treat the adia- 
batic potential self-consistently; later calculations with 
approximate adiabatic functionals suggest having self- 
consistency decreases the error somewhat. 

In Fig. 5, we plot snapshots of the density at 5au in- 
tervals until T ~ 20au for the two choices of exact initial 
states discussed above, evolving in the fixed adiabati- 
cally exact potential plotted in Fig. 4. As each time we 
show the exact static excited-state density. While both 
choices show density Teakrng' out, the $* choice is more 
successful in preventing this, as might be expected from 
the above discussions, but cannot stop the density melt- 
ing away in its outer regions. The $^^[71*] choice is poor 
from the start, jettisoning density outwards and becom- 
ing too narrow. 

This error in the adiabatically exact evolution is 
caused entirely by initial-state dependence: for the anal- 
ogous calculation beginning in the interacting ground- 
state and starting with the ground-state KS wavefunc- 
tion, the adiabatic approximation would be exact at all 
times when there are no perturbing fields. 




FIG. 6: The exact adiabatic xc (solid) and AEXX (dashed) po- 
tentials evaluated on the true excited-state density n* . (Note 
that the solid line here is the same as the dotted line in the top 
panel of Fig. 4 
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FIG. 7: The ground and first excited orbitals from ground-state 
EXX and LDA, compared to those from the exact GS KS poten- 
tial. Also shown are those orbitals found in Sec. IV A of this 
configuration that yield the exact excited-state density. 



C. Approximate GS functionals 

Finally we look at the third source of error, namely us- 
ing an approximate ground-state xc functional instead 
of the exact one. This error contributes at both the 
ground-state level in obtaining the initial orbitals, and 
in the adiabatic functional when computing dynamics. 
We first turn our attention to the latter and investigate 
the adiabatic exact-exchange functional, AEXX, which, 
for a two-electron system, is simply v^^^^[n]{r,t) = 
—Vh [n] (r, t)/2. In Fig. 6 we compare the adiabatically ex- 
act xc potential found in the previous section to AEXX, 
at the initial time. The undulations of the exact adiabatic 
xc potential that are missing in the AEXX represent cor- 
relation; they are relatively small in this particular case 
once added to the external potential, so we expect that 
AEXX dynamics fairly approximates the adiabatically 
exact dynamics in this situation at least for short times. 

We now return to the error of not having the correct 
initial density (Sec. IV A). 

In a usual TDDFT calculation, a ground-state DFT cal- 
culation is first performed and the orbitals from this are 
used to create the initial wavefunction. So we would 
start with the ground-state KS density of Fig. 2 except 



here there is a further error that will be made, as an ap- 
proximate ground-state xc functional is used instead of 
the exact one. 

In Fig. 7 we plot the ground- and first-excited- KS 
orbitals from ground-state DFT calculations using EXX 
and LDA for the soft-Coulomb helium atom, and com- 
pare to the exact KS orbitals. We also plot, for com- 
pleteness, the pair of orbitals found for this configu- 
ration that yield the exact interacting excited-state den- 
sity, i.e. these are the exact orbitals which, when singly- 
occupied, make up the densities shown in Fig. 2. 

If we first compare the approximate orbitals to the ex- 
act KS orbitals of the soft-Coulomb potential, we find 
both LDA and EXX perform very well, particularly for 
the ground-state orbital. Their ground-state KS poten- 
tials are however quite different: it is well known that 
the LDA potential is much too shallow, with a much too 
rapid exponential decay away from the nucleus, while 
the EXX is much closer to the exact, with the correct 
asymptotic decay. This is reflected in the first-excited 
orbital in LDA, which is more diffuse than EXX and the 
first-excited orbital of the exact soft-Coulomb KS po- 
tential. The asymptotic behavior will be important es- 
pecially when we turn on an electric field, so we will 
(mostly) use the EXX orbitals to build our initial KS 
wavefunctions. 

In the next section we will add the spin-symmetry- 
broken initial state (c) introduced in Section II to our 
investigations; while no longer a spin eigenstate like 
the exact case or $* and $ss^ i\^[^ gj-jjj ]-,^g ^-j^g same to- 
tal density as the exact system. We will simply con- 
struct this from the same EXX orbitals composing $*, 
but will evolve the two orbitals using different xc po- 
tentials, f|c and Uxc- For exact-exchange, we have 
i'HxcN(a^) = ^;H[?^^^](a;) for two electrons. We are 
however not doing spin-DFT since the individual spin- 
densities are not physical ones; we only ever consider 
their sum as observable. (Indeed, even the initial spin- 
densities are wrong). The point of considering such an 
initial state and dynamics is that, for two electrons, it 
is an example of dynamics in orbital-specific potentials, 
as would occur in generalized KS approaches [36-38]. 
Orbital-dependent functionals require an optimized ef- 
fective potential approach to find a single potential for 
all the KS orbitals to live in. This procedure is nu- 
merically very intensive, so often a generalized KS ap- 
proach is used, that relaxes the condition that all orbitals 
evolve under the same potential; moreover, the latter 
has shown to have certain advantages over the OEP 
in certain cases, e.g. better band-gaps. The interest in 
orbital-dependent functionals is that that they can work 
less hard to capture both spatially-non-local and time 
non-local- density-dependence since the orbitals them- 
selves are non-local functionals of the density. 

We note nevertheless that it can be shown there is no 
ISD for two-electron dynamics with spin-TDDFT; this 
follows straightforwardly generalizing the arguments of 
Ref. 27 to each spin-density This is not the case for our 
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FIG. 8: The exact dipole moment and those from TDDFT cal- 
culations starting in the self -consistent ground-state KS wave- 
function and propagated with the corresponding adiabatic ap- 
proximation, all under the influence of the electric field de- 
scribed in the text. 




FIG. 9: The dipole moments when beginning in $* composed 
of EXX orbitals (dashed line) and the exact orbitals of Sec. IV A 
(solid line) in the electric field (see text). 



spin-broken case, however, as here we are not attempt- 
ing to reproduce the exact spin-densities but rather the 
total density evolution; only n{r) = rif + ni will be con- 
sidered meaningful. 

Having delineated and explored the three possible 
sources of error in an adiabatic TDDFT calculation of 
excited state dynamics, we are now ready to run a typi- 
cal adiabatic TDDFT calculation on our model systems, 
starting with our different choices of initial KS states, 
and compare with the exact result. 



D. Propagation in an electric field 

We simulate electric-field driven dynamics in the soft- 
Coulomb helium model, using adiabatic TDDFT. We ap- 
ply a relatively weak oscillating electric field field of am- 
plitude O.Olau with an off-resonant frequency of 0.2a.u. 
We run for 6-cycles including a trapezoidal envelope 
consisting of a 2-cycle linear switch on and a 2-cycle 
switch off. 

As a point of reference, we first look at the perfor- 
mance of TDDFT when starting in the ground state. We 









l'^\ 




1 \ 






- 














r 


/ _\ 




1 '■/■' 


a 


4 


'/^ 


v.. 


U 


- 


— EXACT 






^J 














■■■ *SB 
















' 


- 


._.. ^_« 














'^J 


- 









100 



FIG. 10: The dipole moments for dynam^ics in the given electric 
field starting in the $*, $sb, ^^'^ initial KS wavefunctions and 
using AEXX, compared to exact. 
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FIG. 11: Dynamics for the spin-broken case, $sb/ for differ- 
ent orbital-specific adiabatic functionals, EXX, LSD ASIC, and 
LSD A, compared to the exact. 



propagate the doubly-occupied EXX ground-state or- 
bital in the electric field, using AEXX, and compare to 
the exact dynamics of the initial interacting ground- 
state. The resulting dipole moment is very accurate as 
can be seen in Fig. 8. In fact, almost exact d5n-iam- 
ics can be achieved in this case with an adiabatic ap- 
proximation: the same figure shows the result of using 
the adiabatic self -interaction corrected LDA approxima- 
tion (LDASIC) to propagate the LDASIC ground-state 
orbital. The dipole moment lies practically on top of the 
exact curve. Including correlation in a self-interaction 
free functional therefore improves over exact exchange, 
but an adiabatic approximation is certainly adequate in 
this case. We also show for comparison the LDA result. 

The adiabatic approximation is however not so rosy, 
as we shall see, in the case of an initially excited state. 
The initial interacting state is the first singlet excited 
state {E = -1.705au). 

The typical TDDFT calculation will begin in an ex- 
cited state determined by the orbitals of the corre- 
sponding ground-state KS system, but as discussed in 
Sec. IV A, this does not have quite the right density to 
start with. Since in our model system, we found an ini- 
tial KS state of the correct density, we first check the 



error that using the AEXX orbitals make. We plot in 
Fig. 9 the dipole moments for AEXX calculations start- 
ing in the $* initial state with exact orbitals and with 
the ground-state EXX orbitals (i.e. using the solid and 
dotted orbitals of Figure 7 respectively). As 

previously anticipated, the two calculations perform 
similarly, especially at shorter times. 

At last we are ready to study the excited-state elec- 
tron dynamics, using, as in practical calculations, an ap- 
proximate ground-state KS potential (in this case EXX) 
to generate the orbitals composing the various initial 
state choices. In Fig. 10, we plot the dipole mo- 
ment for each choice of initial state. The spin-broken- 
symmetry case gives the best dynamics; as mentioned 
earlier, functionals expressed directly in terms of instan- 
taneous orbitals automatically contain time non-local 
and spatially-non-local density-information. In this par- 
ticular case the orbital-specific EXX potentials, obtained 
from spin-scaling the (spin-restricted) AEXX potential, 
are self-interaction-free, while the latter is not for the 
case of the excited two-orbital singlet state. 

Unlike in the propagation of the ground-state, adding 
correlation to the adiabatic potential does not signifi- 
cantly improve the results, as can be seen in Fig. 11; 
the symmetry-broken LSDASIC and EXX calculations 
are both off by roughly the same amount. Fig. 11 also 
shows the LSDA result, clearly worse than the others by 
comparison. 

The $* initial state has the correct spin symmetry and 
is closest in character to the true excited state, being a 
double Slater determinant. It gives better dynamics than 
the $s''[n*] initial state, whose dynamics are completely 
incorrect. This is consistent with the results in the previ- 
ous sections. We can say with confidence that it is ISD, 
and its underlying static correlation that it causing poor 
dynamics for this doubly-occupied-singlet case. 

We can understand the poor result for the initial state 
of the Eq. (8) form also by looking at the exact adiabatic 
xc potential and Vxc[n*,^*,^^^[n*]]{x), both shown in 
Fig. 4. This form treats the density as if it were a ground- 
state and uses one doubly occupied orbital, inversion of 



this orbital yields 
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n\{x), i.e. this is the potential 



in which the wavefunction is static at the initial time. 
However in a calculation using the adiabatic approxi- 
mation , it will fall into a much deeper potential, dras- 
tically effecting its dynamics. Without the initial-state 
dependence of the xc potential giving a bump similar to 
Fig. 1, it will perform poorly. 



V. MODEL LIH DYNAMICS 

Our final example is perhaps the most topical one 
where the initial state is an excited one: coupled 
electron-ion d5mamics after photo-excitation. The ini- 
tial excitation is assumed to place the initial nuclear 
wavepacket on an excited potential energy surface, ver- 
tically up from its ground-state equilibrium. Field-free 




FIG. 12: The exact Born-Oppenheimer potential energy sur- 
faces (solid lines) and the EXX KS Born-Oppenheimer surfaces 
(dashed lines). Note that due to degeneracy between the low- 
est two surfaces, these calculations cannot be converged for 
large values of R, the interatomic separation. 
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FIG. 13: The interatomic distance, R{t), for Ehrenfest dynam- 
ics starting in the first excited-state and propagated freely. 
The various KS initial wavefunctions ($*, ^sb) ^^^) are prop- 
agated with the adiabatic exact exchange approximation. In- 
sert: short-time dynamics for the various initial states. 



dynamics on the excited surface ensues; usually the nu- 
clei are treated classically, coupled in either an Ehren- 
fest or surface-hopping scheme to the quantum elec- 
tron dynamics. For Ehrenfest dynamics in the TDDFT 
framework [33], and the simplest type of surface- 
hopping [39], the electrons start in the KS excited-state 
obtained simply from promoting an electron from the 
highest occupied orbital in the ground-state KS config- 
uration to a virtual orbital. This initial state is then 
evolved in the time-dependent external potential caused 
by the moving ions. The dynamics of the ions is given by 
simple Newtonian mechanics where the electronic sys- 
tem provides a force / n(r, i)VHWoxt(r; {R-})/ where {R} 
represent the nuclear coordinates. All the three errors 
mentioned earlier therefore raise their heads: the initial 
density is not that of the true excited state (even if the ex- 
act xc potential was used), an adiabatic approximation 
is used to propagate it, and this involves an approxi- 
mate ground-state functional. It should be noted that in 
the more accurate surface-hopping method of Ref. [40], 
where correct TDDFT surfaces are used to provide the 
nuclear forces [41], the first problem does not arise; the 
method in fact circumvents having to define the elec- 
tronic state. 

We will use the parameters found in Ref. 42 for a 
one-dimensional two-electron lithium hydride model, 
where 



Vc^x) = -l/v/a;2 + 0.7- 1/Va;2 + 2.25 (18) 



and Vnn{R) = -1.0/\/i?2 + 1.95. We integrate Newton's 
equations of motion using the leapfrog algorithm, and 
for simplicity, we use the same time step for both elec- 
trons and ions. 

The Born-Oppenheimer potential energy surfaces are 
show^n in Fig. 12 for both the exact system and those 
of the bare KS system calculated within spin restricted 
EXX. 

The ground state surface has a minimum at i? = 1.55, 
whereas the EXX has the gs equilibrium at i? = 1.45. 



Near equilibrium the surfaces are somewhat similar. 
The spin-restricted EXX surfaces however encounter the 
well-known fractional charge problem as R increases, 
with the calculations becoming extremely difficult to 
converge and the lowest two surfaces collapsing onto 
each other. In the exact surface, an avoided crossing 
can be seen between the ground and first excited state 
around R = A, where the latter establishes ionic charac- 
ter, eventually decaying as —0.896 — 1/x, while the neu- 
tral dissociation curve flattens out. In fact there are mul- 
tiple avoided crossings at larger distances, the next be- 
ing around i? = 10. Unrestricted EXX calculations break 
the spin-symmetry of the true ground-state at a critical 
separation, but describe dissociation and the shape of 
the potential energy surfaces at larger distances qualita- 
tively [43, 44] . 

We vertically excite the system at the ground-state 
equilibrium geometry into the first excited singlet state 
and then begin the propagation with the initial nuclear 
momentum chosen as zero. In Fig. 13, we show the in- 
ternuclear separation as a function of time for our three 
choices of initial states, calculated within EXX. For $* 
we make the usual practise of taking i/io in Eqs. 7 and 
10 as the HOMO and 0i as the LUMO; likewise for the 
spin-symmetry-broken version '^gg. The density of this 
state n* is then used to define the doubly-occupied or- 
bital in the SSD state <i>«-^[n*] (Eq. 8). 

For this particular system, the exact system remains 
on the first excited surface; the nuclear coordinate os- 
cillates between the turning points at i? = 1.5 and 
R = 11.165. The ground state ^^^ completely fails to 
capture this behavior, making only very small oscilla- 
tions around its initial position. The excited states <f>* 
and <I>gg show reasonable agreement for short times, 
with the spin-broken state once again edging out the 
other. For longer times however, the spin-broken state is 
dramatically wrong, yielding nuclei moving away from 
each other with constant speed. The excited state <f>* 
shows behavior more qualitatively exact but, because of 
the deviation of its potential energy curve compared to 
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the exact one, the oscillation period is quite wrong. 

The conclusions regarding ISD in this example, at 
least at short times, are consistent with the results 
throughout this paper and support the idea that when 
using an adiabatic approximation choosing a KS initial- 
state with a configuration close to that of the exact in- 
teracting initial state yields the better result. However, 
the behavior at longer times in this problem highlights, 
above all, the need for accurate ground-state function- 
als: EXX produces poor nuclear dynamics when $* is 
chosen as the KS state primarily because the shape of 
its potential energy surfaces is wrong - the ground state 
as well as the excited state ones. Any adiabatic approx- 
imation that does not have strong correlation will not 
perform well. Despite giving improved potential energy 
surfaces, the spin-broken EXX approach fails for dynam- 
ics for longer times. We believe this is because as the 
molecule dissociates, symmetry-breaking localizes each 
electron on one nucleus or the other, with each evolv- 
ing according to a different xc potential; each atom sees 
just a neutral atom, leading to a net zero potential, and 
a constant velocity (somewhat as if the molecule was 
dissociating on the ground-state surface). It should be 
noted that the spin-broken solution is not evolving on 
the spin-broken PES. e.g. even at the initial time, at the 
equilibrium geometry there is no symmetry-breaking in 
a spin-unrestricted approach, so the spin-broken surface 
is on top of the spin-restricted surface showed in the 
figure. However the spin-broken wavefunction (c) is 
evolved using different potentials for each orbital, un- 
like the evolution dictated by either the spin-restricted 
EXX surface or the unrestricted EXX (the symmetry- 
breaking point occurs around R = 2.5). 



VI. CONCLUSIONS AND OUTLOOK 

Until relatively recently, ISD fell in the realm of a the- 
oretical curiosity, but due to an increasing number of 
topical applications beginning with the system not in its 
ground-state, ISD is now also of much practical concern. 
Almost all calculations today use adiabatic functionals 
that neglect the ISD that the exact functionals are known 
to have. 

By considering several choices of initial KS states 
when the initial interacting state is excited, we explored 
the effect of ISD on dynamics in several exactly-solvable 
model systems and the performance of the adiabatic ap- 
proximation. We noted there are three sources of error: 
excited KS initial states do not have the same density 
as the corresponding true excited states even when the 
exact functional is used, the use of the adiabatic approx- 



imation to propagate the initial state, and the ground- 
state functional generating the adiabatic approximation 
itself being approximate. By separating these errors 
from each other, we were able to properly assign how 
badly they impacted the dynamics in model systems, al- 
lowing us to see the effect of ISD more clearly. 

When the initial KS state has a configuration that is 
significantly different than the true state, the adiabatic 
approximation fails severely - even for non-interacting 
electrons, where there is a significant error that can be 
interpreted as a static correlation effect. The optimal 
choice for the KS initial state is one whose configuration 
is most similar to the true excited state; this is in fact 
the usual practise in recent calculations (e.g. Ref. [33]). 
The error in using an adiabatic xc functional to propa- 
gate such a state, however, is significantly larger than 
the error the same functional produces when describ- 
ing the propagation of an initial true ground-state, as 
demonstrated by the model soft-Coulomb helium atom 
example. For two electrons, the spin-broken initial 
KS state, with the two orbitals evolved under different 
spin-decomposed potentials (depending on the instan- 
taneous spin-densities), often appeared to be the best 
choice, suggesting that orbital-specific potentials (al- 
lowed by the generalized KS scheme) could be a useful 
future avenue of research for TDDFT. However, how to 
obtain such potentials in a general iV-electron case is not 
obvious (at least without the need for empirical param- 
eters). Generally, orbital-dependent functionals treated 
within the OEP may be promising for these problems 
where memory-(of the density)-dependence, including 
ISD at the KS level only, is naturally captured by the 
KS orbitals; again, devising suitable orbital-dependent 
functionals with the appropriate ISD is a direction for 
future research. 

The example of coupled electron-nuclear dynamics 
in the model LiH molecule, although supporting the 
earlier results of the optimal choice of a KS initial 
state, above all, highlighted the need for more accu- 
rate ground-state functional approximations for disso- 
ciation. 

We do not provide in this work an approximation that 
includes ISD, both of the true interacting state and the 
KS state; this is likely a difficult task. But what is clear 
from the results in this paper, and the stage of the field, 
is that now is the right time to address this issue. 
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